Source code for hysop.backend.device.opencl.operator.spatial_filtering

# Copyright (c) HySoP 2011-2024
#
# This file is part of HySoP software.
# See "https://particle_methods.gricad-pages.univ-grenoble-alpes.fr/hysop-doc/"
# for further info.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.


import numpy as np
import functools

from hysop.tools.htypes import check_instance, first_not_None
from hysop.tools.decorators import debug
from hysop.tools.numpywrappers import npw
from hysop.backend.device.opencl.opencl_operator import OpenClOperator
from hysop.core.graph.graph import op_apply
from hysop.fields.continuous_field import Field
from hysop.parameters.parameter import Parameter
from hysop.topology.cartesian_descriptor import CartesianTopologyDescriptors
from hysop.operator.base.spatial_filtering import (
    RemeshRestrictionFilterBase,
    SpectralRestrictionFilterBase,
    SubgridRestrictionFilterBase,
    PolynomialInterpolationFilterBase,
    PolynomialRestrictionFilterBase,
)
from hysop.backend.device.opencl.opencl_symbolic import OpenClSymbolic
from hysop.backend.device.opencl.opencl_copy_kernel_launchers import (
    OpenClCopyBufferRectLauncher,
)
from hysop.backend.device.opencl.opencl_kernel_launcher import OpenClKernelListLauncher
from hysop.backend.device.opencl.opencl_elementwise import (
    OpenClElementwiseKernelGenerator,
)
from hysop.symbolic import local_indices_symbols
from hysop.symbolic.relational import Assignment


[docs] class OpenClPolynomialInterpolationFilter( PolynomialInterpolationFilterBase, OpenClOperator ):
[docs] @debug def discretize(self): if self.discretized: return super().discretize() dFin = self.dFin dFout = self.dFout gr = self.grid_ratio dim = dFin.dim assert dFin.is_scalar assert dFout.is_scalar assert self.subgrid_interpolator.gr == gr ekg = self.elementwise_kernel_generator Wr = self.subgrid_interpolator.Wr n = self.subgrid_interpolator.n ghosts = np.asarray(self.subgrid_interpolator.ghosts) I = np.asarray(local_indices_symbols[:dim][::-1]) fin, fout = ekg.dfields_to_ndbuffers(dFin, dFout) Fin = ekg.symbolic_tmp_scalars("F", shape=n, dtype=dFin.dtype) Fout_values = Wr.dot(Fin.ravel()).reshape(gr) exprs = () for idx in np.ndindex(*n): e = Assignment(Fin[idx], fin(I + idx - ghosts)) exprs += (e,) for idx in np.ndindex(*gr): e = Assignment(fout(gr * I + idx), Fout_values[idx]) exprs += (e,) kname = f"interpolate_grid_{self.polynomial_interpolation_method}".lower() interpolate_grid_kernel, _ = ekg.elementwise_kernel( kname, *exprs, compute_resolution=self.iter_shape, debug=False ) exchange_ghosts = self.dFout.exchange_ghosts(build_launcher=True) kl = OpenClKernelListLauncher(name=kname, profiler=self._profiler) kl += interpolate_grid_kernel kl += exchange_ghosts self.execute_kernels = functools.partial(kl, queue=self.cl_env.default_queue)
@op_apply def apply(self, **kwds): super().apply(**kwds) evt = self.execute_kernels()
[docs] class OpenClPolynomialRestrictionFilter( PolynomialRestrictionFilterBase, OpenClOperator ):
[docs] @debug def discretize(self): if self.discretized: return super().discretize() dFin = self.dFin dFout = self.dFout gr = self.grid_ratio dim = dFin.dim assert dFin.is_scalar assert dFout.is_scalar assert self.subgrid_restrictor.gr == gr ekg = self.elementwise_kernel_generator Rr = self.subgrid_restrictor.Rr / self.subgrid_restrictor.gr ghosts = np.asarray(self.subgrid_restrictor.ghosts) I = np.asarray(local_indices_symbols[:dim][::-1]) fin, fout = ekg.dfields_to_ndbuffers(dFin, dFout) def gen_inputs(*idx): return fin(gr * I + idx - ghosts) input_values = np.asarray( tuple(map(gen_inputs, np.ndindex(*Rr.shape))) ).reshape(Rr.shape) output_value = (Rr * input_values).sum() e = Assignment(fout(I), output_value) exprs = (e,) kname = f"restrict_grid_{self.polynomial_interpolation_method}".lower() restriction_grid_kernel, _ = ekg.elementwise_kernel( kname, *exprs, compute_resolution=self.iter_shape, debug=False ) exchange_ghosts = self.dFout.exchange_ghosts(build_launcher=True) kl = OpenClKernelListLauncher(name=kname, profiler=self._profiler) kl += restriction_grid_kernel kl += exchange_ghosts self.execute_kernels = functools.partial(kl, queue=self.cl_env.default_queue)
@op_apply def apply(self, **kwds): super().apply(**kwds) evt = self.execute_kernels()
[docs] class OpenClSubgridRestrictionFilter(SubgridRestrictionFilterBase, OpenClSymbolic): """ OpenCL implementation for lowpass spatial filtering: small grid -> coarse grid using the subgrid method. """ def __init__(self, **kwds): super().__init__(**kwds) Fin = self.Fin Fout = self.Fout dim = Fin.dim assert Fin.is_scalar assert Fout.is_scalar # We do not know the grid ratio and array strides before discretization. # so we defer the initialization of those integers with symbolic constants. (symbolic_input_buffer,) = self.symbolic_buffers("fine_grid") symbolic_output_buffer = self.Fout.s() symbolic_grid_ratio = self.symbolic_constants("gr", count=dim, dtype=npw.int32) symbolic_input_strides = self.symbolic_constants( "is", count=dim, dtype=npw.int32 ) symbolic_input_ghosts = self.symbolic_constants( "gs", count=dim, dtype=npw.int32 ) I = local_indices_symbols[:dim][::-1] read_idx = npw.dot( symbolic_input_strides, npw.add(npw.multiply(symbolic_grid_ratio, I), symbolic_input_ghosts), ) expr = Assignment(symbolic_output_buffer, symbolic_input_buffer[read_idx]) self.require_symbolic_kernel("extract_subgrid", expr) self.symbolic_input_buffer = symbolic_input_buffer self.symbolic_output_buffer = symbolic_output_buffer self.symbolic_grid_ratio = symbolic_grid_ratio self.symbolic_input_strides = symbolic_input_strides self.symbolic_input_ghosts = symbolic_input_ghosts
[docs] @debug def setup(self, work): dFin, dFout = self.dFin, self.dFout ibuffer, obuffer = dFin.sbuffer, dFout.sbuffer self.symbolic_input_buffer.bind_memory_object(ibuffer) for i in range(dFin.dim): self.symbolic_grid_ratio[i].bind_value(self.grid_ratio[i]) self.symbolic_input_strides[i].bind_value( ibuffer.strides[i] // ibuffer.dtype.itemsize ) self.symbolic_input_ghosts[i].bind_value(dFin.ghosts[i]) super().setup(work) (extract_subgrid, _) = self.symbolic_kernels["extract_subgrid"] exchange_ghosts = self.dFout.exchange_ghosts(build_launcher=True) kl = OpenClKernelListLauncher(name="extract_subgrid") kl += extract_subgrid kl += exchange_ghosts self.execute_kernels = functools.partial(kl, queue=self.cl_env.default_queue)
@op_apply def apply(self, **kwds): evt = self.execute_kernels()
[docs] class OpenClSpectralRestrictionFilter(SpectralRestrictionFilterBase, OpenClOperator): """ OpenCL implementation for lowpass spatial filtering: small grid -> coarse grid using the spectral method. """ def _compute_scaling_coefficient(self): kernel_generator = OpenClElementwiseKernelGenerator( cl_env=self.cl_env, kernel_config=self.kernel_config ) # Kernels to copy src_slices to dst_slices (windowing operation on frequencies) kl = OpenClKernelListLauncher(name="lowpass_filter") for src_slc, dst_slc in zip(*self.fslices): kl += OpenClCopyBufferRectLauncher.from_slices( "copy", src=self.FIN[src_slc], dst=self.FOUT[dst_slc] ) # Now we compute the scaling coefficient of the filter # self.Ft.input_buffer is just a pypencl.Array so we need to use # the kernel_generator to fill ones and use explicit copy kernels # This seems to be the only solution to fill a non C-contiguous # OpenClinput_buffer with ones. buf = self.Ft.input_buffer (buf,) = kernel_generator.arrays_to_symbols(buf) expr = Assignment(buf, 1) fill_ones, _ = kernel_generator.elementwise_kernel("fill_ones", expr) fill_ones(queue=self.cl_env.default_queue) self.Ft(simulation=False) kl(queue=self.cl_env.default_queue) # here we apply unscaled filter self.Bt(simulation=False) # Here we get the coefficient scaling = 1.0 / self.Bt.output_buffer[(0,) * self.FOUT.ndim].get() # Now we can finally build the filter scaling kernel (fout,) = kernel_generator.arrays_to_symbols(self.FOUT) expr = Assignment(fout, scaling * fout) scale, _ = kernel_generator.elementwise_kernel("scale", expr) kl += scale # we store the filtering kernel list for the setup step self.filter = functools.partial(kl, queue=self.cl_env.default_queue) # finally build ghost exchanger exchange_ghosts = self.dFout.exchange_ghosts(build_launcher=True) if exchange_ghosts is not None: exchange_ghosts = functools.partial( exchange_ghosts, queue=self.cl_env.default_queue ) self.exchange_ghosts = exchange_ghosts return scaling @op_apply def apply(self, **kwds): """Apply spectral filter (which is just a square window centered on low frequencies).""" super().apply(**kwds) evt = self.Ft(**kwds) evt = self.filter() evt = self.Bt(**kwds) if self.exchange_ghosts is not None: evt = self.exchange_ghosts()